Zamień sygnał na sumę sygnałów (1 pkt) np
x = sin.(2*pi*t*200) + 2* sin.(2*pi*t*400)
Zaobserwuj wynik transformaty i wyjaśnij go.
using Plots
using FFTW
function plot_pressure(time_array, signal, title)
return plot(time_array, signal, title=title, xlabel="Time [ms]", ylabel="Pressure", label="")
end
function plot_fourier_transform(signal_fft, title)
return sticks(abs.(signal_fft), title=title, xlabel="Frequency [Hz]", ylabel="Power", label="")
end
function scale_wave_to_rms(signal, required_rms)
rms = sqrt(mean(signal.^2))
scale_factor = required_rms / rms
return signal * scale_factor
end;
# sampling frequency
Fs = 1024
# time points in interval [0, 1] every 1/Fs
t = 0:1/(Fs-1):1;
x = sin.(2*pi*t*200) + 2 * sin.(2*pi*t*400);
plot_pressure(t, x, "Signal")
y = fft(x);
plot_fourier_transform(y, "Fourier transform")
sticks(real(y))
sticks(imag(y))
The real-input (r2c) DFT in FFTW computes the forward transform $Y$ of the size $n$ real array $X$, exactly as defined above, i.e. $$Y_k=\sum_{j=0}^{n-1}X_je^{-2\pi ijk/n}$$ This output array $Y$ can easily be shown to possess the “Hermitian” symmetry $Y_k = Y_{n-k}^*$, where we take $Y$ to be periodic so that $Y_n = Y_0$.
As a result of this symmetry, half of the output $Y$ is redundant (being the complex conjugate of the other half), and so the 1d r2c transforms only output elements $0…\frac{n}{2}$ of $Y$ ($\frac{n}{2+1}$ complex numbers), where the division by 2 is rounded down.
This property arises because the complex exponential term $e^{-2\pi ijk/n}$ in the DFT equation has the property $e^{-2\pi ij(n-k)/n} = e^{-2\pi ijk/n}*$, which leads to the Hermitian symmetry.
Usuwanie szumów (1 pkt) :
1. Wypełniamy tablicę wartościami funkcji cosinus ("sygnału") zaburzonej niewielkim "szumem" np. dodając do każdej wartości wylosowaną liczbę funkcją rand().
2. Proszę narysować wykres zaszumionej funkcji.
3. Narysować wykres transformaty Fouriera (widmo) tego sygnału (jak poprzednio).
4. Po transformacie wyzerowac w widmie wszystkie elementy, których wartość bezwzględna jest mniejsza niz 50. W ten sposób usuwamy "szumy" z sygnału.
5.Przeprowadzić odwrotną transformatę funkcją ifft(). Narysować wykres otrzymanej funkcji dla <b>częsci rzeczywistej</b>. Porównać z wejściowym wykresem sygnału.
# sampling frequency
Fs = 1024
# time points in interval [0, 1] every 1/Fs
t = 0:1/(Fs-1):1;
x = cos.(2*pi*t*200)+rand(Float64, length(t));
plot_pressure(t, x, "Signal with noise")
y = fft(x);
plot_fourier_transform(y, "Fourier transform of signal with random noise")
corrected_y = [abs(z) < 50 ? 0.0+0.0im : z for z in y];
typeof(corrected_y)
Vector{ComplexF64} (alias for Array{Complex{Float64}, 1})
plot_fourier_transform(corrected_y, "Corrected Fourier transform")
corrected_x = ifft(corrected_y);
plot_pressure(t, real(corrected_x), "Corrected signal")
W poprawionym sygnale wyraźnie widać początkową funkcje cosinus, szczególnie w środkowym przedziale.
Proszę nagrać własny glos i zastosować na nim trasformatę Fouriera, narysować wykres widma. Następnie poeksperymentować (wyciąć wybrane częstotliwości), dokonać odwrotnej transformaty i odsłuchać efekt (3 pkt) .
Przydatne materiały: 1. Basic sound processing 2. Pakiet Wav
using Pkg
#Pkg.add("WAV")
using WAV
using Plots
using Statistics
# Used for making some sound louder
# sound, sampling_frequency = wavread("my_voice.wav");
# size(sound)
# sound *= 8
# wavwrite(sound, "my_voice.wav", Fs=sampling_frequency)
# wavplay("my_voice.wav")
# loading requested sound
sound, sampling_frequency = wavread("my_voice.wav");
n, _ = size(sound)
(13312, 2)
# we'll work only with chanel 1 from now onwards
s1 = sound[:,1];
wavwrite(s1, "Original.wav", Fs=sampling_frequency)
wavplay("Original.wav")
# Time of recording [s]
sound_time = size(s1, 1) / sampling_frequency
0.6037188f0
# time array
time_array = (0:(size(s1, 1)-1)) / sampling_frequency
time_array *= 1000 # scaling to miliseconds
0.0f0:0.04535147f0:603.67346f0
plot_pressure(time_array, s1, "Original sound")
# Root mean square - measure of waveform amplitude
rms_val = sqrt(mean(s1.^2))
0.04147038109257873
s1_fft = fft(s1);
plot_fourier_transform(s1_fft, "Original sound")
title = "Low frequencies removed";
f_bound = 2000
s1_fft_modified = copy(s1_fft)
s1_fft_modified[1:f_bound] = zeros(ComplexF64, f_bound);
s1_fft_modified[n-f_bound+1:n] = zeros(ComplexF64, f_bound);
plot_fourier_transform(s1_fft_modified, title)
s1_modified = real(ifft(s1_fft_modified))
# Scale the modified waveform back to original power
s1_modified = scale_wave_to_rms(s1_modified, rms_val);
plot_pressure(time_array, s1_modified, title)
wavwrite(s1_modified, title*".wav", Fs=sampling_frequency)
wavplay("Original.wav")
wavplay(title*".wav")
title="High frequencies removed";
f_bound = 1000
s1_fft_modified = copy(s1_fft)
s1_fft_modified[f_bound:n-f_bound-1] = zeros(ComplexF64, n-2*f_bound);
plot_fourier_transform(s1_fft_modified, title)
s1_modified = real(ifft(s1_fft_modified))
# Scale the modified waveform back to original power
s1_modified = scale_wave_to_rms(s1_modified, rms_val);
plot_pressure(time_array, s1_modified, title)
wavwrite(s1_modified, title*".wav", Fs=sampling_frequency)
wavplay("Original.wav")
wavplay(title*".wav")
title="Low power frequencies removed";
f_bound = 5
s1_fft_modified = copy(s1_fft)
s1_fft_modified = [abs(z) < f_bound ? 0.0+0.0im : z for z in s1_fft_modified];
plot_fourier_transform(s1_fft_modified, title)
s1_modified = real(ifft(s1_fft_modified))
# Scale the modified waveform back to original power
s1_modified = scale_wave_to_rms(s1_modified, rms_val);
plot_pressure(time_array, s1_modified, title)
wavwrite(s1_modified, title*".wav", Fs=sampling_frequency)
wavplay("Original.wav")
wavplay(title*".wav")
title="High power frequencies removed";
f_bound = 15
s1_fft_modified = copy(s1_fft)
s1_fft_modified = [abs(z) > f_bound ? 0.0+0.0im : z for z in s1_fft_modified];
plot_fourier_transform(s1_fft_modified, title)
s1_modified = real(ifft(s1_fft_modified))
# Scale the modified waveform back to original power
s1_modified = scale_wave_to_rms(s1_modified, rms_val);
plot_pressure(time_array, s1_modified, title)
wavwrite(s1_modified, title*".wav", Fs=sampling_frequency)
wavplay("Original.wav")
wavplay(title*".wav")
title="Stronger low frequencies";
addend = 1
multiplier = 50
n
13312
We will smoothly, in linear manner, multiply low frequencies from 1 up to f_bound
scale_function = x -> (abs(x/n - 1/2)) * multiplier + addend
aux = [i for i in 1:n];
plot(aux, scale_function.(aux), title="Scale function", label="")
s1_fft_modified = copy(s1_fft)
s1_fft_modified = [z*scale_function(i) for (i, z) in enumerate(s1_fft_modified)];
plot_fourier_transform(s1_fft_modified, title)
s1_modified = real(ifft(s1_fft_modified))
# Scale the modified waveform back to original power
s1_modified = scale_wave_to_rms(s1_modified, rms_val);
plot_pressure(time_array, s1_modified, title)
wavwrite(s1_modified, title*".wav", Fs=sampling_frequency)
wavplay("Original.wav")
wavplay(title*".wav")
title="Stronger high frequencies";
addend = 0.5
multiplier = 100
n
13312
scale_function = x -> (-(abs(x/n - 1/2)) + 1/2) * multiplier + addend
aux = [i for i in 1:n];
plot(aux, scale_function.(aux), title="Scale function", label="")
s1_fft_modified = copy(s1_fft)
s1_fft_modified = [z*scale_function(i) for (i, z) in enumerate(s1_fft_modified)];
plot_fourier_transform(s1_fft_modified, title)
s1_modified = real(ifft(s1_fft_modified))
# Scale the modified waveform back to original power
s1_modified = scale_wave_to_rms(s1_modified, rms_val);
plot_pressure(time_array, s1_modified, title)
wavwrite(s1_modified, title*".wav", Fs=sampling_frequency)
wavplay("Original.wav")
wavplay(title*".wav")